##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmpi2uNTb/downloaded_packages
##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmpi2uNTb/downloaded_packages
##
## The downloaded binary packages are in
## /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmpi2uNTb/downloaded_packages
options(future.globals.maxSize = 74 * 1024^3) # 55 GB
getOption("future.globals.maxSize") #59055800320## [1] 79456894976
load(here("jk_code", "JK_cleanMD.rds"))
SO4 <- subset(SO4, features = grep("^mt-|^Rp|^Gm", rownames(SO4), invert = TRUE, value = TRUE))
SO4 <- SCTransform(SO4) %>%
RunPCA() %>%
FindNeighbors(dims = 1:15) %>%
FindClusters(resolution = 0.1) %>%
RunUMAP(dims = 1:15)## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 15822 by 11431
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
## Found 286 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 15822 genes
## Computing corrected count matrix for 15822 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 16.47543 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Set default assay to SCT
## PC_ 1
## Positive: Umod, Egf, Tmem52b, Sult1d1, Sostdc1, Fabp3, Klk1, Prdx5, Foxq1, Cox6c
## Wfdc15b, Ly6a, Wfdc2, Krt7, Slc25a5, Ckb, Cox5b, Atp5g1, Cldn19, Cox4i1
## Ggt1, Atp1a1, Ndufa4, Cox8a, Atp5b, Chchd10, Gadd45g, Cox6b1, Atp1b1, Cox7b
## Negative: Pappa2, Zfand5, Aard, Robo2, Pde10a, Wwc2, Itga4, Nadk2, Nos1, S100g
## Ramp3, Neat1, Sgms2, Ptgs2, Col4a4, Irx1, Col4a3, Tmem158, Mir6236, Ranbp3l
## Itprid2, Bmp3, Cdkn1c, Camk2d, Sdc4, Mcub, Dctd, Etnk1, Srrm2, Peg3
## PC_ 2
## Positive: Fos, Junb, Jun, Hspa1b, Hspa1a, Btg2, Egr1, Zfp36, Atf3, Ier2
## Fosb, Socs3, Klf6, Klf2, Jund, Dnajb1, Dusp1, Tsc22d1, Gadd45g, Rhob
## H3f3b, Gadd45b, Cebpd, Ubc, Mt1, Mt2, Actb, H2bc4, H1f2, Klf4
## Negative: Mir6236, Pappa2, CT010467.1, Egf, Slc12a1, Umod, Etnk1, Wnk1, Sfrp1, Atp1b1
## Nme7, Robo2, Sptbn1, Col4a4, Hsp90b1, Sec14l1, Pde10a, Zbtb20, Itprid2, Dst
## Mal, App, Lars2, Wwc2, Col4a3, Kcnq1ot1, Atrx, Utrn, Tfcp2l1, Pou3f3
## PC_ 3
## Positive: Fth1, Ubb, Ftl1, Ldhb, Car15, Fxyd2, Prdx1, Cd63, Cd9, Eif1
## Mpc2, Mgst1, Mt1, Clu, Bsg, Tmem213, Aldoa, Mdh1, Ramp3, Itm2b
## Spp1, S100a1, Selenow, Tmem59, Tmem176b, Wfdc2, Tmbim6, Tspo, Atpif1, Ppp1r1a
## Negative: Mir6236, Egf, CT010467.1, Umod, Neat1, Tmem52b, Nme7, Fos, Malat1, Slc12a1
## Kcnq1ot1, Jun, Junb, Lars2, Etnk1, Dst, Egr1, Wnk1, Slc5a3, Sult1d1
## Foxq1, Atrx, Atp1b1, Pnisr, Fosb, Atf3, Ivns1abp, Zfp36, Btg2, Syne2
## PC_ 4
## Positive: Pappa2, Aard, Tmem52b, Umod, Egf, Sult1d1, Tmem158, Foxq1, Mcub, Ptgs2
## Dctd, Wwc2, Iyd, Car15, Ramp3, Ptprz1, Hsp90b1, Cd9, S100g, Wnk1
## Cdkn1c, Defb42, 1700028P14Rik, Pth1r, 5330417C22Rik, Tmsb4x, Tagln2, Ctsc, Clu, Bmp2
## Negative: Mt1, Apoe, Mt2, Aebp1, Sostdc1, Gpx6, Fxyd2, Ptger3, Neat1, Ckb
## Fgf9, Egfl6, Ivns1abp, Mfsd4a, Defb1, Car4, Plet1, CT010467.1, Fkbp11, Chchd10
## Atp5md, 2900052N01Rik, Cox6c, Igfbp5, Atp5k, Tmem213, Chka, Mgst3, Avpr1a, Atp1a1
## PC_ 5
## Positive: Hspa1a, Hspa1b, CT010467.1, Pappa2, Klk1, Car15, Cldn10, Fth1, Hspa8, Jun
## Lars2, Fau, Mir6236, Itm2b, Wfdc2, Uba52-ps, Aard, Ftl1, Ptger3, Eef1a1
## Pik3r1, Egf, Sfrp1, Mal, Gpx4, Tspo, Atp1a1, Tmem176a, Id3, Hsp90aa1
## Negative: S100g, Actb, Tmem52b, Sdc4, Tmsb10, Ppia, Abhd2, Uroc1, Egr1, Serf2
## Slc39a1, Igfbp7, Atf3, Ndufb1-ps, Atp5md, Alkbh5, Cebpd, Cox6c, Ramp3, Ndufa3
## Gnb1, Atp5e, Foxq1, Atp5k, Ldhb-ps, Fam107a, Kdm6b, Rbm47, Atp5mpl, Wfdc15b
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11431
## Number of edges: 355780
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9249
## Number of communities: 4
## Elapsed time: 1 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 17:42:46 UMAP embedding parameters a = 0.9922 b = 1.112
## 17:42:46 Read 11431 rows and found 15 numeric columns
## 17:42:46 Using Annoy for neighbor search, n_neighbors = 30
## 17:42:46 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 17:42:46 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//Rtmpi2uNTb/file11f4f16975343
## 17:42:46 Searching Annoy index using 1 thread, search_k = 3000
## 17:42:48 Annoy recall = 100%
## 17:42:48 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 17:42:49 Initializing from normalized Laplacian + noise (using RSpectra)
## 17:42:49 Commencing optimization for 200 epochs, with 480882 positive edges
## 17:42:49 Using rng type: pcg
## 17:42:51 Optimization finished
SO4@meta.data <- SO4@meta.data %>%
mutate(subclass_MD = dplyr::case_when(
seurat_clusters == 0 ~ "type_1",
seurat_clusters == 1 ~ "type_2",
seurat_clusters == 2 ~ "type_3",
seurat_clusters == 3 ~ "type_4"
))
SO4@meta.data$subclass_MD <- factor(SO4@meta.data$subclass_MD , levels = c("type_1", "type_2", "type_3", "type_4"))
Idents(SO4) <- SO4@meta.data$subclass_MD
DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)## Calculating cluster type_1
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
## Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster type_2
## Calculating cluster type_3
## Calculating cluster type_4
SO.markers %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
slice_head(n = 10) %>%
ungroup() -> top10
DoHeatmap(SO4, features = top10$gene) + NoLegend()## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Ifi47,
## Rsad2, Gbp9, Gbp3, Oasl2
type_1_markers <- (SO.markers[SO.markers$cluster == "type_1", ])
type_2_markers <- (SO.markers[SO.markers$cluster == "type_2", ])
type_3_markers <- (SO.markers[SO.markers$cluster == "type_3", ])
type_4_markers <- (SO.markers[SO.markers$cluster == "type_4", ])#Type 1
df<- type_1_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>% rownames_to_column(var="SYMBOL")
markers <- DEG_list %>%
rownames_to_column(var = "SYMBOL")
ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 1.68% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC))) %>%
arrange(p_val_adj) %>%
head(100)
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:100] "22682" "268902" "56089" "68646" "18125" "23984" "16371" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...329 enriched terms found
## 'data.frame': 329 obs. of 12 variables:
## $ ID : chr "GO:0001822" "GO:0072001" "GO:0072073" "GO:0051926" ...
## $ Description : chr "kidney development" "renal system development" "kidney epithelium development" "negative regulation of calcium ion transport" ...
## $ GeneRatio : chr "13/97" "13/97" "8/97" "6/97" ...
## $ BgRatio : chr "375/28928" "390/28928" "173/28928" "76/28928" ...
## $ RichFactor : num 0.0347 0.0333 0.0462 0.0789 0.0519 ...
## $ FoldEnrichment: num 10.34 9.94 13.79 23.54 15.46 ...
## $ zScore : num 10.56 10.31 9.79 11.41 9.77 ...
## $ pvalue : num 4.15e-10 6.68e-10 1.28e-07 2.20e-07 3.73e-07 ...
## $ p.adjust : num 8.11e-07 8.11e-07 1.03e-04 1.33e-04 1.67e-04 ...
## $ qvalue : num 5.57e-07 5.57e-07 7.10e-05 9.16e-05 1.15e-04 ...
## $ geneID : chr "Robo2/Irx1/Epcam/Col4a4/Col4a3/Irx2/Sdc4/Cdkn1c/Bcl2/Hoxb7/Enpp1/Sfrp1/Nrp1" "Robo2/Irx1/Epcam/Col4a4/Col4a3/Irx2/Sdc4/Cdkn1c/Bcl2/Hoxb7/Enpp1/Sfrp1/Nrp1" "Robo2/Irx1/Epcam/Irx2/Sdc4/Bcl2/Hoxb7/Sfrp1" "Nos1/Mcub/Bcl2/Calm1/Ptgs2/Calm2" ...
## $ Count : int 13 13 8 6 7 9 7 9 6 6 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 lines## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
df<- type_2_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>%
rownames_to_column(var = "SYMBOL")
ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 2.95% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC))) %>%
arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:50] "330428" "16513" "53315" "268379" "15220" "69824" "13645" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...97 enriched terms found
## 'data.frame': 97 obs. of 12 variables:
## $ ID : chr "GO:0006119" "GO:0009060" "GO:0045333" "GO:0019646" ...
## $ Description : chr "oxidative phosphorylation" "aerobic respiration" "cellular respiration" "aerobic electron transport chain" ...
## $ GeneRatio : chr "16/50" "17/50" "17/50" "12/50" ...
## $ BgRatio : chr "154/28928" "206/28928" "271/28928" "66/28928" ...
## $ RichFactor : num 0.1039 0.0825 0.0627 0.1818 0.1481 ...
## $ FoldEnrichment: num 60.1 47.7 36.3 105.2 85.7 ...
## $ zScore : num 30.6 28 24.3 35.3 31.8 ...
## $ pvalue : num 7.88e-25 1.27e-24 1.48e-22 7.82e-22 1.10e-20 ...
## $ p.adjust : num 6.18e-22 6.18e-22 4.80e-20 1.90e-19 2.14e-18 ...
## $ qvalue : num 4.32e-22 4.32e-22 3.35e-20 1.33e-19 1.50e-18 ...
## $ geneID : chr "Cox5b/Uqcrq/Cox6c/Cox5a/Chchd10/Uqcr10/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Ndufb9/Cox4i1/Cox7c/Uqcrh/Cox8a/Ndufa13" "Idh2/Cox5b/Uqcrq/Cox6c/Cox5a/Chchd10/Uqcr10/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Ndufb9/Cox4i1/Cox7c/Uqcrh/Cox8a/Ndufa13" "Idh2/Cox5b/Uqcrq/Cox6c/Cox5a/Chchd10/Uqcr10/Uqcr11/Ndufb8/Cox6b1/Cox7a2/Ndufb9/Cox4i1/Cox7c/Uqcrh/Cox8a/Ndufa13" "Cox5b/Uqcrq/Cox5a/Uqcr10/Uqcr11/Ndufb8/Cox7a2/Ndufb9/Cox4i1/Cox7c/Uqcrh/Cox8a" ...
## $ Count : int 16 17 17 12 12 12 17 12 12 6 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 li## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
df<- type_3_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>%
rownames_to_column(var = "SYMBOL")
ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 24.08% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC))) %>%
arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:50] "14282" "13653" "12702" "11910" "14281" "16598" "16477" "22695" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...451 enriched terms found
## 'data.frame': 451 obs. of 12 variables:
## $ ID : chr "GO:1903706" "GO:0009408" "GO:0035914" "GO:0042026" ...
## $ Description : chr "regulation of hemopoiesis" "response to heat" "skeletal muscle cell differentiation" "protein refolding" ...
## $ GeneRatio : chr "11/48" "6/48" "6/48" "4/48" ...
## $ BgRatio : chr "458/28928" "96/28928" "105/28928" "24/28928" ...
## $ RichFactor : num 0.024 0.0625 0.0571 0.1667 0.0317 ...
## $ FoldEnrichment: num 14.5 37.7 34.4 100.4 19.1 ...
## $ zScore : num 11.9 14.7 14 19.9 11 ...
## $ pvalue : num 1.85e-10 1.25e-08 2.15e-08 6.92e-08 7.79e-08 ...
## $ p.adjust : num 3.00e-07 1.01e-05 1.16e-05 2.10e-05 2.10e-05 ...
## $ qvalue : num 1.54e-07 5.20e-06 5.95e-06 1.08e-05 1.08e-05 ...
## $ geneID : chr "Fos/Zfp36/Jun/Hspa1b/Tsc22d1/Zfp36l1/Snai2/Hsp90aa1/Irf1/Nfkbiz/Runx1" "Hspa1a/Hspa1b/Dnajb1/Hsp90aa1/Ier5/Dnaja1" "Egr1/Atf3/Fos/Btg2/Egr2/Nr4a1" "Hspa1a/Hspa1b/Hsp90aa1/Dnaja1" ...
## $ Count : int 11 6 6 4 7 7 8 7 4 4 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 li## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
df<- type_4_markers %>% arrange(desc(avg_log2FC))
df2 <- df %>% filter(p_val_adj < 0.05)
DEG_list <- df2
markers <- DEG_list %>%
rownames_to_column(var = "SYMBOL")
ENTREZ_list <- bitr(
geneID = rownames(DEG_list),
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db
)## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 37.5% of input gene IDs are fail to map...
markers <- ENTREZ_list %>% inner_join(markers, by = "SYMBOL")
# Removing genes that are not statistically significant.
markers <- markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)
pos.markers <- markers %>% dplyr::filter(avg_log2FC > 0) %>% arrange(desc(abs(avg_log2FC))) %>%
arrange(p_val_adj) %>% head(50)
#head(pos.markers, n = 50)
pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)
pos_go <- enrichGO(gene = pos.ranks, #a vector of entrez gene id
OrgDb = "org.Mm.eg.db",
ont = "BP",
readable = TRUE) #whether mapping gene ID to gene Name
pos_go## #
## # over-representation test
## #
## #...@organism Mus musculus
## #...@ontology BP
## #...@keytype ENTREZID
## #...@gene chr [1:25] "15957" "23962" "55932" "626578" "236573" "58185" "100039796" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...104 enriched terms found
## 'data.frame': 104 obs. of 12 variables:
## $ ID : chr "GO:0035456" "GO:0035458" "GO:0009615" "GO:0042832" ...
## $ Description : chr "response to interferon-beta" "cellular response to interferon-beta" "response to virus" "defense response to protozoan" ...
## $ GeneRatio : chr "12/24" "11/24" "13/24" "8/24" ...
## $ BgRatio : chr "77/28928" "68/28928" "412/28928" "43/28928" ...
## $ RichFactor : num 0.1558 0.1618 0.0316 0.186 0.0365 ...
## $ FoldEnrichment: num 188 195 38 224 44 ...
## $ zScore : num 47.3 46.1 21.8 42.2 22.6 ...
## $ pvalue : num 1.36e-25 1.26e-23 1.78e-18 8.63e-18 9.17e-18 ...
## $ p.adjust : num 6.52e-23 3.03e-21 2.85e-16 8.82e-16 8.82e-16 ...
## $ qvalue : num 3.49e-23 1.62e-21 1.53e-16 4.73e-16 4.73e-16 ...
## $ geneID : chr "Ifit1/Gbp3/Tgtp2/Ifi47/Iigp1/Tgtp1/Ifit3/Gbp2/F830016B08Rik/Gbp7/Ifitm3/Irgm1" "Ifit1/Gbp3/Tgtp2/Ifi47/Iigp1/Tgtp1/Ifit3/Gbp2/F830016B08Rik/Gbp7/Irgm1" "Ifit1/Oasl2/Rsad2/Tgtp1/Nlrc5/Oasl1/Rtp4/Ifit3b/Ifit3/Gbp2/Gbp7/Ifitm3/Irgm1" "Gbp3/Gbp10/Gbp9/Tgtp2/Iigp1/Tgtp1/Gbp2/Gbp7" ...
## $ Count : int 12 11 13 8 12 10 8 8 6 4 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
dotplot(pos_go) +
ggtitle("") +
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "left",
axis.text.y = element_text(hjust = 0, size = 10)) +
scale_y_discrete(position = "right",
labels = function(x) str_wrap(x, width = 25)) # Wrap y-axis labels to 2 li## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
gene_list <- c("Ptgs2", "Nos1", "Mcub", "Egf",
"Cox5b",
"Fos", "Egr1", "Jun",
"Cxcl10","Ifit1")
DoHeatmap(SO4, features = gene_list) + NoLegend()